# -*- coding: utf-8 -*-
"""
Created on Wed Apr 10 14:13:35 2013
TODO: change class of performance and mission results, provide easy in/out
@author: Maxim
"""
import aircraft
import aerodynamics
import FlightConditions as FC
#import performance
from performance import *
from numpy import array, linspace, hstack, zeros, log, sin,cos,tan, arcsin, ones,cumsum
from scipy.integrate import simps
from miscTools import Timer
import matplotlib.pyplot as plt
import engine
timer = Timer()
from propeller_design import FixedPitchThrust

class MissionResults:
    """
    stores results of mission analysis such as velocity, mass, power setting etc. 
    at different segments of the mission.
    """
    def __init__(self,nseg=10):
        self.nseg = nseg
        self.n    = 0
        self.totalMass       = zeros(nseg)
        self.fuelMass        = zeros(nseg)
        self.velocity        = zeros(nseg)
        self.EAS             = zeros(nseg)
        self.density         = zeros(nseg)
        self.densityAltitude = zeros(nseg)
        self.fuelFlow        = zeros(nseg)
        self.rpm             = zeros(nseg)
        self.powerSetting    = zeros(nseg)
        self.thrust          = zeros(nseg)
        self.power           = zeros(nseg)
        self.propEfficiency  = zeros(nseg)
        self.beta            = zeros(nseg)
        self.SAR_km_per_kg   = zeros(nseg)
        self.drag            = zeros(nseg)
        self.lift            = zeros(nseg)
        self.range           = 0.0
        self.endurance       = 0.0
        self.fuelBurned      = 0.0
    def add_performance_result(self,rslt,fuelMass,totalMass):
        i = self.n
        self.n += 1
        self.velocity[i]        = rslt.velocity
        self.EAS[i]             = rslt.equivalentAirspeed
        self.densityAltitude[i] = rslt.densityAltitude
        self.density[i]         = rslt.density
        self.fuelFlow[i]        = rslt.fuelFlow
        self.rpm[i]             = rslt.RPM
        self.powerSetting[i]    = rslt.powerSetting
        self.thrust[i]          = rslt.thrust
        self.power[i]           = rslt.power
        self.propEfficiency[i]  = rslt.propEfficiency
        self.beta[i]            = rslt.beta
        self.SAR_km_per_kg[i]   = rslt.SAR_km_per_kg
        self.drag[i]            = rslt.drag
        self.lift[i]            = rslt.lift
        self.fuelMass[i]        = fuelMass
        self.totalMass[i]       = totalMass
        if self.nseg>=2:
            self._integrate()
    def _integrate(self):
        self.range = simps(self.SAR_km_per_kg[:self.n],self.fuelMass[:self.n])
        specificEndurance = 1.0/(array(self.fuelFlow[:self.n])*3600.0)
        self.endurance = simps(specificEndurance,self.fuelMass[:self.n])
        self.fuelBurned = abs(self.fuelMass[0] - self.fuelMass[self.n])
    def _disp(self,name,value,unit=''):
        return '{0:18}{1:8.2f} {2:6}'.format(name,value,unit)
    def _disp_sequence(self,name,value,unit=''):
        value = value[:self.n]
        out = '{0:18}'.format(name)
        for val in value:
            out += '{0:8.2f}'.format(val)
        return out + ' {0:6}'.format(unit)
    def display(self):
        title = 'Mission results'
        print title + '\n' + '-'*len(title)
        print self._disp('Range',self.range,'km')
        print self._disp('Endurance',self.endurance,'h')
        print self._disp('Fuel burned',self.fuelBurned,'kg')
        print self._disp_sequence('Power setting',self.powerSetting)
        print self._disp_sequence('Beta setting',self.beta,'deg')
        print self._disp_sequence('rpm',self.rpm)
        print self._disp_sequence('fuel mass',self.fuelMass,'kg')


class Cruise(SteadyLevelFlight):
    def get_reserve_fuel(self,time=30.0,powerSetting=75.0):
        """
        Calculates mass of reserve fuel
        
        Parameters
        ----------
        time : float, min
            time to stay in air
        """
        rslt = MissionResults()
        time = time*60.0
        rslt.endurance = time
        self.run_velocity_maxTAS(self.basicInput, self.thrustModule,powerSetting)
        fuelMass = self.fuelFlow*time
        rslt.fuelBurned = fuelMass
        rslt.range = self.velocity*time
        print '{0:.2f} kg for {1:.1f} min flight at {2:.0f}% power setting'.format(fuelMass,time/60.0,powerSetting)
        return rslt
    def get_maximum_range_breguet(self,startFuelmass,endFuelmass,altitude,n=3):
        """
        Calculates maximum range using Breguet equation
        """
        extInp = ExtendedInput(self.ac,2.1,altitude)
        extInp.update_fuel_mass(startFuelmass)
        w0 = extInp.ac.mass()
        extInp.update_fuel_mass(endFuelmass)
        w1 = extInp.ac.mass()
        extInp.update_fuel_mass((startFuelmass+endFuelmass)/2.0)
        self.run_velocity_maxSAR(extInp,self.thrustModule)
        rslt = MissionResults()
        effy = self.propEfficiency
        c = self.fuelFlow
        LD = self.L_D
        rslt.range = effy/c * LD * log(w1/w0)
        return rslt
    def get_maximum_range(self,startFuelmass,endFuelmass,altitude,n=3):
        """
        Calculates maximum range using integration wrt fuel weight
        """
        rslt = MissionResults()
        fuelWeight = linspace(endFuelmass,startFuelmass,n)
        for w in fuelWeight:
            self.bi.update_fuel_mass(w)
            self.run_velocity_maxSAR()
            totalMass = self.ac.mass()
            rslt.add_performance_result(self,w,totalMass)
        return rslt
    def get_range_at_power_setting(self,startFuelmass,endFuelmass,altitude,powerSet,n=3):
        """
        Calculates range at given power setting.
        """
        rslt = MissionResults()
        fuelWeight = linspace(endFuelmass,startFuelmass,n)
        for w in fuelWeight:
            self.bi.update_fuel_mass(w)
            self.run_velocity_maxTAS(powerSet)
            rslt.add_performance_result(self,w,self.bi.ac.mass())
        return rslt



class Climb(ClimbAndDescendingFlight):
    def _init_climb(self,aircraft,basicInput,thrustModule):
        self.ac = aircraft
        self.bi = basicInput
        self.tm = thrustModule
    def climb_simulation_max(self,startAltitude,endAltitude,nStep=5,powerSetting=95.0):
        # climb simulation using direct methods
        altitude = linspace(startAltitude,endAltitude,nStep)
        fuelMass = 0.0
        timeToClimb = 0.0
        distance = 0.0
        dh = altitude[2] - altitude[1]
        for h in altitude[1:]:
            timer.start()
            self.bi.atm.set_altitude(h)
            self.bi.density = self.bi.atm.density
            self.run_maximumClimbRate(self.bi,self.tm, powerSetting)
            dt = dh / self.climbRate
            dWf = dt*self.fuelFlow
            distance += dt*self.velocity*cos(self.climbAngle)
            fuelMass += dWf
            timeToClimb += dt
            timer.finish()
        print distance, fuelMass, timeToClimb
        return distance, fuelMass, timeToClimb
    def climb_simulation_integr(self,startAltitude,endAltitude,nStep=5,powerSetting=95.0):
        # climb simulation using numerical integration
        # one of the methods should be selected for future development
        altitude = linspace(startAltitude,endAltitude,nStep)
        dt, ff, V = zeros(nStep),zeros(nStep),zeros(nStep)
        for i,h in enumerate(altitude):
            timer.start()
            self.bi.atm.set_altitude(h)
            self.bi.density = self.bi.atm.density
            self.run_maximumClimbRate(self.bi,self.tm,powerSetting)
            if i==0:
                dt[i] = 0.0
            else:
                dt[i] = (altitude[i] - altitude[i-1])/self.climbRate
            ff[i] = self.fuelFlow
            V[i] = self.velocity
            timer.finish()
        distance = simps(V,dt)
        fuelBurned = simps(ff,dt)
        timeToClimb = sum(dt)
        return distance,fuelBurned,timeToClimb
            
    def glide_simulation(self):
        pass

class Takeoff:
    """
    Takeoff analysis class. Takeoff consists of 3 segments: ground roll, 
    rotation and climbout. All these segments are analyzed separately. 
    Result of analysis is distance, time and fuel burned for each stage of 
    takeoff.
    
    Parameters
    ----------
    aircraft : aircraft object
    altitude : float, m
    """
    def __init__(self,aircraft,altitude,thrustModule,CLto=0.7,CLmaxTo=1.9,
                 hObstacle=15.0):
        self.CLto = CLto
        self.CLmaxTo = CLmaxTo
        self.hObstacle = hObstacle
        self.bi = ExtendedInput(aircraft,CLmaxTo,altitude)
        self.fm = flightMechanics()
        self.tm = aircraft.analysis.thrust
        self.ac = aircraft
        self.alt = altitude
        self.atm = FC.ISAtmosphere(altitude)
        self.rho = self.atm.density
        self._get_data()
        self.mu = 0.2
        self.distance = zeros(3)
        self.time = zeros(3)
        self.fuelWeight = zeros(3)
    
    def analyze(self,dtRot=1.5,powerSet=100.0):
        self.distance[0], self.time[0], self.fuelWeight[0] = self._ground_roll(powerSet)
        self.distance[1], self.time[1], self.fuelWeight[1] = self._rotation(powerSet,dtRot)
        self.distance[2], self.time[2], self.fuelWeight[2] = self._climbout(powerSet,dtRot)
        self.totalDistance = sum(self.distance)
        self.totalTime = sum(self.time)
        self.totalFuelWeight = sum(self.fuelWeight)
        return self.totalDistance, self.totalTime, self.totalFuelWeight

    def _get_data(self):
        hw = self.ac.wing.aapex[2] - self.ac.landingGear.groundContact_Z[2]
        b = self.ac.wing.span
        hw16 = (16.0*hw/b)**2.0
        self.Ge = hw16 / (1.0+hw16) # ground effect
        self.Vstall = self.fm.vStall(self.bi)
        self.Vlo = 1.1 * self.Vstall
        self.Vtr = 1.2 * self.Vstall
        self.Vflare = (self.Vlo + self.Vtr) / 2.0

    def _ground_roll(self,powerSet,nSeg=10):
        tol = 1e-5
        err = tol+1.0
        g = self.bi.g
        CDto = self.bi.Cd0 + self.Ge*self.bi.k*self.CLto**2.0
        V = linspace(0.0,self.Vlo,nSeg)
        dV = V[1]-V[0]
        Q = 0.5 * self.bi.density * V**2.0
        L = Q*self.bi.refArea*self.CLto
        D = Q*self.bi.refArea*CDto
        T = zeros(nSeg)
        ff= zeros(nSeg)
        dt= ones(nSeg)
        for i,v in enumerate(V):
            self.tm.runAnalysis(v,self.bi.density,powerSet)
            T[i] = self.tm.thrust
            ff[i] = self.tm.fuelFlow
        while err>tol:
            tLast = sum(dt)
            Wf = g*dV*cumsum(ff*dt)
            W = self.bi.wt*ones(nSeg)-Wf
            F = T-D-self.mu*(W-L)
            dt = (W/g)/F
            err = abs(tLast-sum(dt))
        ds = V*dt
        t = cumsum(dt*dV)
        s = cumsum(ds*dV)
        Sgr = s[-1]
        tgr = t[-1]
        Wfgr = W[0]-W[-1]
        return Sgr, tgr, Wfgr
        
    def _rotation(self,powerSet,dtRot=1.5):
        Srot = self.Vlo*dtRot
        self.tm.runAnalysis(self.Vlo,self.bi.density,powerSet)
        wf = self.bi.g * self.tm.fuelFlow*dtRot
        return Srot, dtRot, wf

    def _climbout(self,powerSet,dtRot=1.5):
        Wtr = self.bi.wt - sum(self.fuelWeight)
        nFlare = self.bi.density*self.Vflare**2/2.0 * self.bi.refArea*self.CLmaxTo / Wtr
        R = self.Vflare / nFlare
        self.tm.runAnalysis(self.Vtr,self.bi.density,powerSet)
        fuelFlow = self.tm.fuelFlow
        #Wfrot = self.bi.g*fuelFlow*dtRot
        Q = self.bi.density*self.Vtr*self.Vtr / 2.0
        CL = Wtr / (Q*self.bi.refArea)
        CD = self.bi.get_drag_coefficient(CL)
        D = Q*self.bi.refArea*CD
        climbAngle = arcsin((self.tm.thrust-D)/Wtr)
        #-->
        Sflare = R*climbAngle
        hFlare = R*(1.0-cos(climbAngle))
        tFlare = Sflare / self.Vflare
        WfuelFlare = fuelFlow*tFlare
        #-->
        hTR = self.hObstacle - hFlare
        RCtr = self.Vtr*sin(climbAngle)
        ttr = hTR/RCtr
        Str = self.Vtr*cos(climbAngle)*ttr
        WfTR = self.bi.g*fuelFlow*ttr
        
        Sclimb = Sflare + Str
        timeClimb = tFlare + ttr
        WfClimb = WfuelFlare + WfTR
        return Sclimb, timeClimb, WfClimb


class Landing():
    def __init__(self,aircraft,altitude,thrustModule,CLgr=0.7,CLmaxLDG=2.2,
                 hObstacle=15.0,mu=0.40):
        self.ac = aircraft
        self.bi = ExtendedInput(aircraft,CLmaxLDG,altitude)
        
        self.hObst = hObstacle
        self.atm = FC.ISAtmosphere(altitude)
        self.fm = flightMechanics()
        self.tm = thrustModule
        self.alt = altitude
        self.rho = self.atm.density
        self.reactionTime = 1.0
        self.mu = mu
        self.idlePower,self.idleFuelFlow = aircraft.engine.engineList[0].get_idle_parameters()
        self.idlePower = 0.0
        self.loadFactor = 1.2
        self.CLgr = CLgr
        self.CLmaxLdg = CLmaxLDG
        self._get_data()
        self.distance = zeros(2)
        self.time = zeros(2)
        self.fuelWeight = zeros(2)
    def analyze(self):
        self.distance[0],self.time[0],self.fuelWeight[0] = self._approach()
        self.Wgr = self.bi.wt  -self.fuelWeight[0]*self.bi.g
        self.distance[1],self.time[1],self.fuelWeight[1] = self._ground_roll()
        return sum(self.distance),sum(self.time),sum(self.fuelWeight)

    def _get_data(self):
        hw = self.ac.wing.aapex[2] - self.ac.landingGear.groundContact_Z[2]
        b = self.ac.wing.span
        hw16 = (16.0*hw/b)**2.0
        self.Ge = hw16 / (1.0+hw16) # ground effect
        self.Vstall = self.fm.vStall(self.bi)
        self.Vap = 1.3*self.Vstall
        self.Vtd = 1.1*self.Vstall
        self.Va = (self.Vap+self.Vtd) / 2.0
        self.W = self.ac.mass() * self.bi.g
    
    def _get_cl(self,V):
        return 2.0*self.W / (self.atm.density*V*V*self.ac.wing.area)
    
    def _approach(self):
        T = self.idlePower * self.Va / 2.0
        ff = self.idleFuelFlow
        CL = self._get_cl(self.Va)
        CD = self.bi.get_drag_coefficient(CL)
        theta = arcsin(1.0/(CL/CD)-T/self.W)
        R = self.Va**2.0 / (self.bi.g*(self.loadFactor - 1.0))
        hf = R*(1.0 - cos(theta))
        Sa = (self.hObst - hf) / tan(theta) + R*sin(theta)
        Vg = self.Va * cos(theta)
        ta = Sa / Vg
        Mfa = ff * ta
        return Sa, ta, Mfa*self.bi.g

    def _ground_roll(self,Nseg=50):
        CDldg = self.bi.Cd0 + self.Ge*self.bi.k * self.CLgr**2.0
        dV = -self.Vtd / (Nseg-1)
        V = linspace(self.Vtd,0.0,Nseg)
        Q = self.bi.atm.density*V*V / 2.0
        L = Q*self.bi.refArea * self.CLgr
        D = Q*self.bi.refArea * CDldg
        tol = 1.0e-6
        err = tol+1.0
        dt = ones(Nseg)
        ff = ones(Nseg)*self.idleFuelFlow
        while err>tol:
            tLast = sum(dt)
            Wf = self.bi.g * cumsum(ff*dt)
            W = self.bi.wt*ones(Nseg) - Wf
            F = -D - self.mu*(W-L) #ma
            a = F * self.bi.g/ W
            dt = dV / a
            err = abs(tLast-sum(dt))
        ds = V*dt
        t = cumsum(dt)
        s = cumsum(ds)
        Sgr = s[-1] + self.reactionTime*self.Vtd
        tgr = t[-1] + self.reactionTime
        Wfgr = Wf[-1]
        Mfgr = Wfgr / self.bi.g
        return Sgr,tgr,Mfgr*self.bi.g
        
# --- test section ---
def print_header(header):
    print '\n'+header
    print '='*len(header)

def run_test1():
    import matplotlib.pyplot as plt
    name = 'V0510-MDO'
    velocity =50.0
    altitude =2000.
    Clmax    =1.4
    ac       =aircraft.load(name)
    rslt = ac.analyze_aero_trim(velocity,altitude)
    ac.analysis.aerodynamics.update(velocity,1.225)
    tm       =ac.analysis.thrust
    bi = BasicInput(ac,Clmax,altitude)

    print_header('Cruise')
    crs = Cruise(bi,tm)
    rslt = crs.get_range_at_power_setting(80.,10.,2000.0,75.0,6)
    rslt.display()
    crs.display()

def run_takeoff():
    name = 'V0510-MDO'
    velocity =50.0
    altitude =2000.
    ac       =aircraft.load(name)
    mass      =ac.mass.totalMass
    rslt = ac.analyze_aero_trim(velocity,altitude)
    ac.analysis.aerodynamics.update(velocity,1.225)
    Cd0      =ac.get_drag(velocity,altitude)
    tm       =ac.analysis.thrust
    
    TO = Takeoff(ac,0.0,tm)
    print TO.analyze(1.5,800.0)

def run_test3():
    ac = aircraft.load('V0510-MDO')
    altitude = 0.0
    dT = 0
    velocity = 50.0
    mass = ac.mass.totalMass
    rslt = ac.analyze_aero_trim(velocity,altitude)
    k        =rslt.k
    Cd0      =ac.get_drag(velocity,altitude)
    S        =ac.wing.area
    Clmax = 1.8
    
    tm       =ac.analysis.thrust
    bi       =basicInput(altitude,dT,Cd0,k,S,mass,Clmax)
    climb = Climb()
    climb._init_climb(ac,bi,tm)
    print climb.climb_simulation_integr(0,2000)
    print climb.climb_simulation_max(0,2000)

if __name__=="__main__":
    run_test1()
    #run_test1()